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ABSTRACT 

We extend the cosmological linear perturbation theory code CMBFAST to closed geometries. 
This completes the implementation of CMBFAST to all types of geometries and allows the user 
to perform an unlimited search in the parameter space of models. This will be specially useful for 
placing confidence limits on cosmological parameters from existing and future data. We discuss 
some of the technical issues regarding the implementation. 
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1. Introduction 

The study of the anisotropies in the Cosmic Microwave Background (CMB) is revolutionizing cosmology. 
Fluctuations across a range of angular scales have been detected and their precision in the measurement of 
the power spectra are improving steadily. Ground based and balloon born experiments are making effort to 
detect the acoustic features in the power spectrum predicted by inflationary models. In the near future the 
MAP satellite and a few years later the Planck satellite will be able to measure the CMB power spectrum 
with sufficient accuracy to constraint several cosmological parameters with a few percent accuracy percent 
accuracy (Jungman et al. 1996| , [Bond et al. 1997| , [Zaldarriaga et al. 1997| ). 



There are two reasons that explain why the CMB is so powerful at constraining the cosmological 
parameters. First is the improvement in the experimental sensitivity that allows the detection of these small 
anisotropies with increasing precision. Second is our ability to accurately calculate the predictions for the 
different world models. These theoretical predictions can be made with high accuracy because the calculation 
of the power spectrum only involves linear perturbation theory. The same fact that makes the anisotropies 
difficult to detect makes them also easy to calculate. 

The most widely used Boltzmann code to calculate anisotropy spectra, CMBFAST, was developed 



by the authors and has become a standard tool freely available to the community Seljak & Zaldarriaga 



1996, Zaldarriaga, Seljak & Bertschinger 1998. The code is based on the integral solution of the transport 
equation for the photons, the Boltzmann equation. The algorithm is much more efficient than the traditional 
way of solving the hierarchy because it takes advantage of the large difference between the distance photons 
have traveled since recombination and the horizon at that time. 

There have been several studies intended to constrain cosmological parameters with the existing CMB 
data |Lineweaver 199^ , [legmark 19991 , [Podelson fc Knox 1999i These studies rely on the speed of CMBFAST 



to search a multidimensional parameter space for the best fitting model and the confidence regions for the 
parameters. They typically evaluate millions of models which would be prohibitely slow with a traditional 
Boltzmann code taking a hundred times longer than CMBFAST to run. So far the available versions of 
CMBFAST did not allow the user to explore models with spatially closed geometries and thus that part of 
parameter space has not been explored with this code (although closed universe Boltzmann codes do exist. 



see e.g. White & Scott 1996). Since the data are pushing us into the corner of parameter space with vanishing 
curvature, any meaningful error determination must also explore the models with positive curvature. It is 
therefore timely to make CMBFAST available also for this class of models. 

A key part of the CMBFAST algorithm is being able to calculate efficiently the radial part of the 
eigenf unctions of the Laplacian, the ultraspherical Bessel functions. In our implementation of open models 
we calculate these functions by integrating their corresponding differential equation. The straightforward 
extension of this algorithm is not possible for closed models. In this paper we explain the technical difficulties 
that arise for closed models and present the solution that is implemented in the new version of CMBFAST 
(V3.0). 



2. The line of sight method and a brief description of geometries 

In spherical coordinates the unperturbed geometry is described by. 



(1) 
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The scale factor is a(T), the coordinate r is conformal time and the radial coordinate x is related to the 
angular radius r by 

r if-i/^sini^i/^x, K>0 
rix) = \ X, K^O (2) 
[ {-K)-^/^smh{-Ky/^X, K<0. 

with K = 7Jg(f2o ~ l)j where flo includes the contribution from both matter and the cosmological constant. 
Note that with our choice of variables we have dr = —dx along radial geodesies. 

To calculate the CMB power spectrum C; we need to evolve the perturbations in the different species 
(CDM, baryons, photons, neutrinos). The perturbations are so small that linear theory suffices. The standard 
method is first to expand all perturbations in eigenvalues of the Laplacian which evolve independently. In 
a flat spacetime this is a Fourier expansion. The wavenumber (3 labels the different eigenmodes. The 
fluctuations produced by each mode can be calculated separately and then all the contributions to the power 
are added, 

Ci - (4^)2 / pdp P(/3)|AT,(/3)r (3) 



We have introduced IS.ti (/3) , the contribution by a mode of wavenumber (3 to the amplitude of multipole 
I of the temperature anisotropies. The primordial power spectrum is denoted by Equation (||) says 

that to obtain the total power in multipole / we add the contributions from all modes labeled by [3. We 
have written this sum as a continuous integral, which applies to flat and open models. For closed spatial 
hypersurfaces the eigenvalues of the Laplacian are discrete so the integral becomes a discrete sum over /3 
with K^^/^P = 3, 4, 5, • • •. The flrst two modes are pure gauge modes and should be excluded from the sum 
Abbott fc Schacfer 198e| . 



In CMBFAST the amplitudes At; are obtained using the integral solution. 



To 



At/= / dT^'^ix)Sif3,r). (4) 
Jo 

The ultraspherical Bessel functions ^'"f^ix) depend on the wavenumber /3 and the radial coordinate along 
the geodesic % = tq — r. The source S{P, r) is a sum of the different terms that can generate temperature 
anisotropies. Among these are gravitational redshift, Doppler shift and intrinsic density fluctuation in the 
photon baryon fluid. The dominant contribution to the source is produced at the last scattering surface. 
The full expression for S together with the equivalent formulas for the CMB polarization can be found in 



Zaldarriaga, Seljak fc Bcrtschinger 199g 



Equation has several advantages over the more traditional way of calculating the anisotropies. The 
main one is that the sampling in /3 needed to compute accurately the source S (essentially the inverse of the 
sound horizon at recombination) is much smaller than that needed to obtain an accurate representation of 
(approximately the inverse of the distance to the last scattering surface). In the traditional Boltzmann 
method one solves directly for At; so one has to sample f3 with the inverse of the distance to the last scattering 
surface. Obtaining the source for each /? involves solving the perturbation equations for the different species 
so a reduction in the number of equations to be solved greatly enhances the performance of the code. 

In essence traditional Boltzmann codes spend most of their time computing Bessel functions. If we had 
a fast method for computing the Bessel functions we could use equation (^) to make a fast algorithm. This 
is what CMBFAST does. In flat geometries the ultraspherical Bessel functions reduce to ordinary spherical 
Bessel functions ji{f3x) which are a function of a single variable x = Px- t^ic implementation of CMBFAST 
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for flat models we use tables with the values of the spherical Bessel functions which are generated in advance 
and stored in disk. 

For open model we are unable to tabulate the ultraspherical Bessel functions in advance because they 
no longer depend on only one variable. They depend separately on (3 and x because the curvature radius of 
the universe introduces another scale that can be used to turn j3 and x into two dimensionless numbers and 
the ultraspherical Bessel functions depend on these two numbers independently. In our implementation of 
CMBFAST for open models we choose to calculate the Bessel function by solving their differential equation, 



,2 + 



rixf 



u' = (5) 



dx^ 

where u^p{x) = r(x)$j3(x)- 

Although at first this may not seem to be the most efficient method (for example, we could have 
used the recurrence relation for the ultraspherical Bessel functions instead), it is actually very efficient for 
our purposes. The point is that we need the function at a succession of neighboring points in time when 
evaluating equation (Q). To get from one point to the next a Runge Kutta integrator only performs a handful 
of operations, as opposed to the recurrence relation calculation which has to start at sufficiently high /. It is 
true that recurrence relation gives the function value at all but since the final CMB spectra are so smooth 
we actually only use every 50th / or so, so this advantage of recurrence relation becomes less significant. 
Thus although the recurrence relation is the most efficient way to obtain a single value of the Bessel function, 
for our purposes solving the differential equation is a better method than the recurrence relation. Another 



method to calculate these Bessel functions based on the WKB approximation has been proposed (Kosowsky 



1999| ). This method is again a fast way to obtain a single value of any since it only involves a few 



operations per <i>^ call. However, for the integral solution which requires many calls at nearby arguments, 
the available implementation of these formulas does not appear to be faster than solving the differential 
equations. 



Closed models 



In Zaldarriaga, Seljak & Bertschinger 1998 we had developed the necessary formalism to solve the 
perturbation equations in both open and closed universes. We derived the Boltzmann hierarchy for both 
temperature and polarization and the integral solution. The other trivial differences between closed and 
open models like the discreteness of the eigenvalues of the Laplacian were also mentioned. The only problem 
not addressed was the calculation of the Bessel functions in closed models. We discuss this issue here. 

Equation (^) is a second order differential equation and thus has two independent solutions. It can be 
viewed as the Schrodinger equation for a particle with energy E = (3^ in a, potential V ^ l{l + l)/r(x)^- 
Thus for X such that Pr{x) > + 1), where the "energy" is greater than the "potential", the solution 
is oscillatory. For /3r(x) < + 1) there is a growing and a decaying solution, oc 7-'+^ and Up oc 
r~', respectively. The growing solution corresponds to ^^{x)- In order to obtain an accurate numerical 
convergence the integration needs to be started in the regime where 3>^(x) is small, which in our case means 
starting the integration close to (3r{x) ~ I- The equation then needs to be evolved in the direction of 
increasing x until recombination, where x ^ '''o- The integration of 'I'^(x) therefore proceeds in the direction 
of decreasing time, opposite to the evolution of Boltzmann, fluid and Einstein equations. If one were to 
start the integration at early time and evolve the Bessel functions to small radial distances x (i-S- present 
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time) the integration would be numerically unstable, as the decaying mode mJj oc r ' would increasingly 
contaminate the solution. 

In closed models the angular distance r(x) is periodic around 9 = K~^^'^x = t^/'^- The potential V{x) is 
a potential well. For parameters such that the maximum value of 9 required in the evaluation of equation (^) 
is greater than 7r/2 the integration of the differential equation (|^) becomes unstable, similarly to integrating 
in the direction of decreasing x in tii6 open model case. As the integration proceeds the solution becomes 
increasingly contamined by the decaying solution which diverges at 9 — t:. In figure we show lines of 
constant in the (r^m, ^^a) plane for z^ax — 1000. The region upwards of the 9 = n /2 line suffers from this 
instability. 

We can use symmetry arguments analogous to those in quantum mechanics to solve the instability 
problem. The potential well V{x) is symmetric around 9 — n/2 and thus the eigenfunctions are either 
symmetric or antisymmetric around this point. Just like in quantum mechanics we find that modes with 
even (odd) (3 — 1 — 1 are symmetric (antisymmetric). In our implementation of CMBFAST for closed models 
we use this fact to obtain the Bessel functions for regions with 9 > 7r/2. We integrate equation (||) up 




Fig. 1. — Contours of constant 9 in the {Um,^A) plane. The line 9 = corresponds to flat models, with 
closed models above that line. In the region between 9 = tt/2 and 9 = n the calculation of the Bessel 
functions becomes unstable and symmetry can be used to determine their values in that regime. The dots 
show the positions of the models we show in figure The two bottom lines also correspond to 6' = 7r/2,7r 
but in an open universe. 



- 6 - 



to 6 = 7r/2 and extend the solution further using the symmetry of the mode. In practice we apply the 
symmetry operation to the source term S{P,t). We construct a symmetric 5"* = {S{(3,t) + S'(/3,t'))/2 and 
an antisymmetric S'" = {S{(3,t) — S{(3,t'))/2 version of the source. Here r and r' are times that have the 
same r{x), related by r' = K^^'^tt — r. We use one of those in (||) according to the parity of the mode being 
calculated and stop the time integration at the = tq — corresponding to 9 = tt/2. We found 

this to be both a satisfactory solution of the instability and an efficient numerical algorithm. 

In figure ||a we show the power spectra obtained using the new CMBFAST version for several closed 
models. We choose models with n^h'^ = 0.237, f^b/i^ = 0.013 and varying = 0.4,0.6,0.8,1.0,1.2. The 
position of the models in the {Qrm^A) plane is shown in figure |l|. The models we show all have the same 
values of ^mh? and flbh^, so the physics of the acoustic peaks is the same for all the models. They only differ 
in the angular diameter distance to the last scattering surface and their amplitude which was normalized to 
COBE. In figure ||b we rescale the I values so that a given physical scale at recombination corresponds to the 
same angular scale. This is achieved using / — > Ir' {xr) /^{x'r)i where r{xR) is the radial comoving distance 
to recombination and primed variables correspond to the values in a fiducial model. We also renormalized 
the amplitude of the spectra so they are no longer COBE normalized. The agreement at high I overall is 
very good, which is a consistency check for the code. 

As we approach 9 = n, r{x) 0, when we are able to see the antipode. As 6' — > tt the Ci spectra is 
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Fig. 2. — Sample of closed universe models calculated with the new version of CMBFAST. Panel b shows 
the same curves but rescaled in I and amplitude to match as closely as possible. 
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compressed so that the acoustic peaks appear at larger and larger scales. When 9 — n all the CMB photons 
that arrive to us from every direction come form one single point, so there are no anisotropics. We have 
restricted the range of models implemented in CMBFAST to those with 9 < tt. 

In our previous implementations of CMBFAST we took advantage of the smoothness of the Ci spectra 
to calculate for every 50th I in the region of the acoustic peaks. This was sufficient to obtain one percent 
accuracy in flat models and became even better for open model were the spectra are stretched further out 
in I. For closed models the situation is reversed because the spectra are compressed so this sampling can 
become insufficient. If the first acoustic peak occurs at / ^ 130(100) as opposed to I ^ 200 for fiat universes, 
the sampling every 50 Is introduces an error of 2% (5%). 



4. Conclusions 

We extended the CMBFAST code to the closed universes. The main technical issue in this implementa- 
tion was the computation of ultraspherical Bessel functions, which was obtained by solving their differential 
equation jut like in the case of open models. We use the symmetry of the ultraspherical Bessel functions 
around 6* = 7r/2 to obtain their values a,t 6 > tt/2 without integration. The code only works for models with 
6 < TT. The new version of CMBFAST for closed universes will allow the exploration of the full parameter 



space. The results of such an exploration will be presented elsewhere Tegmark & Zaldarriaga 1999 
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